In [1]:
%reload_ext autotime
import rasterio
import requests
from rasterio.features import shapes, sieve, rasterize
from rasterio.plot import show
import geopandas as gpd
import os
import numpy as np
from tqdm.auto import tqdm
from shapely.geometry import box, shape
from rasterio.mask import mask
from matplotlib import pyplot as plt
import xdem
import pandas as pd
tqdm.pandas()
pd.set_option('display.max_columns', None)
✔️ 2.09 s (2024-12-02T11:03:52/2024-12-02T11:03:54)
In [2]:
manifest = pd.json_normalize(requests.get("https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/collection.json").json()["links"])
manifest = manifest[manifest["rel"] == "item"]
manifest["tilename"] = manifest.href.str.replace(".json", "").str.strip("./")
manifest
✔️ 211 ms (2024-12-02T11:03:54/2024-12-02T11:03:54)
Out[2]:
| rel | href | type | file:checksum | tilename | |
|---|---|---|---|---|---|
| 2 | item | ./BD43_10000_0204.json | application/json | 12208b52bfdff84d924430e9503bc70bcc7555c1861e5b... | BD43_10000_0204 |
| 3 | item | ./BD43_10000_0205.json | application/json | 12204ec4b1257e965c8da142c4696b1f8cbbc6c78d0cf0... | BD43_10000_0205 |
| 4 | item | ./BD43_10000_0304.json | application/json | 1220acfd4fd86a3866f4facf92bb75cd5224451e774ff4... | BD43_10000_0304 |
| 5 | item | ./BD43_10000_0305.json | application/json | 122051507adf6a042f09c138f1842911fec7bd43c270d8... | BD43_10000_0305 |
| 6 | item | ./BD43_10000_0404.json | application/json | 1220c63b7011af35652df009f8a6b69f4e3f8ea2db33d6... | BD43_10000_0404 |
| ... | ... | ... | ... | ... | ... |
| 302 | item | ./BH43_10000_0102.json | application/json | 12200f2d3649a677e745f258162ca6ebc13904d6ecde95... | BH43_10000_0102 |
| 303 | item | ./BH43_10000_0201.json | application/json | 1220d28f9ff3b2f6ad76139b4ccfa5bd86bc996b417e36... | BH43_10000_0201 |
| 304 | item | ./BH43_10000_0202.json | application/json | 1220d78d97c89b75b5d35667f77acb8e5a2078d52b977b... | BH43_10000_0202 |
| 305 | item | ./BH43_10000_0301.json | application/json | 12200b058bab70e426eccd9ae990a5888e1754b27b97c2... | BH43_10000_0301 |
| 306 | item | ./BH43_10000_0302.json | application/json | 12203c93b413b89941df279575d9d352046cc8b6977deb... | BH43_10000_0302 |
305 rows × 5 columns
In [3]:
tilename = "BD44_10000_0503"
old = rasterio.open(f"https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2018-2020/dem_1m/2193/{tilename}.tiff")
new = rasterio.open(f"https://nz-elevation.s3.ap-southeast-2.amazonaws.com/gisborne/gisborne_2023/dem_1m/2193/{tilename}.tiff")
✔️ 929 ms (2024-12-02T11:03:55/2024-12-02T11:03:55)
In [4]:
new.crs, old.crs
✔️ 2.56 ms (2024-12-02T11:03:56/2024-12-02T11:03:56)
Out[4]:
(CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'),
CRS.from_wkt('PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",GEOGCS["NZGD2000",DATUM["New_Zealand_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6167"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4167"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",173],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",1600000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","2193"]]'))
In [5]:
show(old, title="2018-2020 DEM")
✔️ 4.8 s (2024-12-02T11:03:56/2024-12-02T11:04:01)
Out[5]:
<Axes: title={'center': '2018-2020 DEM'}>
In [6]:
show(new, title="2023 DEM")
✔️ 4.88 s (2024-12-02T11:04:01/2024-12-02T11:04:06)
Out[6]:
<Axes: title={'center': '2023 DEM'}>
In [7]:
diff = new.read(1) - old.read(1)
print(diff.shape)
plt.imshow(diff, cmap="coolwarm_r", vmin=-1, vmax=1)
plt.colorbar()
✔️ 1.4 s (2024-12-02T11:04:06/2024-12-02T11:04:07)
(7200, 4800)
Out[7]:
<matplotlib.colorbar.Colorbar at 0x7f75e42eb430>
In [8]:
_ = plt.hist(diff.flatten(), bins=100, range=(-2, 2))
plt.title("Histogram of elevation differences")
✔️ 448 ms (2024-12-02T11:04:07/2024-12-02T11:04:08)
Out[8]:
Text(0.5, 1.0, 'Histogram of elevation differences')
In [9]:
result = diff.round().clip(min=-1, max=1).astype(np.int16)
result = sieve(result, 4000)
result = np.where(result != 0, result, np.nan)
plt.imshow(result, cmap="coolwarm_r")
plt.colorbar()
✔️ 3.23 s (2024-12-02T11:04:08/2024-12-02T11:04:11)
Out[9]:
<matplotlib.colorbar.Colorbar at 0x7f75e413f820>
In [10]:
areas = gpd.GeoDataFrame(geometry=[shape(s) for s, v in shapes(result, transform=new.transform) if not np.isnan(v)])
areas["area"] = areas.area
areas.sort_values("area", ascending=False, inplace=True)
areas.plot("area", legend=True)
✔️ 1.35 s (2024-12-02T11:04:11/2024-12-02T11:04:12)
Out[10]:
<Axes: >
In [11]:
largest_area = areas.iloc[0]#.head(1)
largest_area
✔️ 6.1 ms (2024-12-02T11:04:12/2024-12-02T11:04:12)
Out[11]:
geometry POLYGON ((2054419 5803910, 2054419 5803909, 20... area 530308.0 Name: 176, dtype: object
In [12]:
masked_old, transform = mask(old, [largest_area.geometry], nodata=np.nan)
nonan = np.argwhere(~np.isnan(masked_old[0]))
top_left = nonan.min(axis=0)
bottom_right = nonan.max(axis=0)
print(top_left, bottom_right)
masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_old)
✔️ 1.25 s (2024-12-02T11:04:13/2024-12-02T11:04:14)
[5290 0] [7044 3057]
Out[12]:
<Axes: >
In [13]:
masked_new, transform = mask(new, [largest_area.geometry], nodata=np.nan)
masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_new)
✔️ 1.02 s (2024-12-02T11:04:14/2024-12-02T11:04:15)
Out[13]:
<Axes: >
In [14]:
features = {
"old_min": np.nanmin(masked_old),
"old_max": np.nanmax(masked_old),
"old_mean": np.nanmean(masked_old),
"old_median": np.nanmedian(masked_old),
"old_std": np.nanstd(masked_old),
"new_min": np.nanmin(masked_new),
"new_max": np.nanmax(masked_new),
"new_mean": np.nanmean(masked_new),
"new_median": np.nanmedian(masked_new),
"new_std": np.nanstd(masked_new),
}
features
✔️ 131 ms (2024-12-02T11:04:15/2024-12-02T11:04:15)
Out[14]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585}
In [15]:
raster = rasterize([largest_area.geometry], out_shape=result.shape, transform=new.transform)
masked_diff = np.where(raster, diff, np.nan)
masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
show(masked_diff)
✔️ 1.82 s (2024-12-02T11:04:15/2024-12-02T11:04:17)
Out[15]:
<Axes: >
In [16]:
features.update({
"diff_min": np.nanmin(masked_diff),
"diff_max": np.nanmax(masked_diff),
"diff_mean": np.nanmean(masked_diff),
"diff_median": np.nanmedian(masked_diff),
"diff_std": np.nanstd(masked_diff),
})
features
✔️ 78.1 ms (2024-12-02T11:04:17/2024-12-02T11:04:17)
Out[16]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585,
'diff_min': -6.8740234,
'diff_max': 43.521973,
'diff_mean': 8.370939,
'diff_median': 4.911972,
'diff_std': 9.06017}
In [17]:
attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
attributes = xdem.terrain.get_terrain_attribute(
masked_diff,
resolution=old.res,
attribute=attribute_names
)
plt.figure(figsize=(8, 6.5))
for i in range(8):
plt.subplot(4, 2, i + 1)
plt.imshow(attributes[i].squeeze(), cmap="viridis")
cbar = plt.colorbar()
cbar.set_label(attribute_names[i])
plt.xticks([])
plt.yticks([])
plt.tight_layout()
plt.show()
✔️ 26.7 s (2024-12-02T11:04:17/2024-12-02T11:04:44)
In [18]:
for i, name in enumerate(attribute_names):
features.update({
f"{name}_min": np.nanmin(attributes[i]),
f"{name}_max": np.nanmax(attributes[i]),
f"{name}_mean": np.nanmean(attributes[i]),
f"{name}_median": np.nanmedian(attributes[i]),
f"{name}_std": np.nanstd(attributes[i]),
})
features
✔️ 655 ms (2024-12-02T11:04:44/2024-12-02T11:04:45)
Out[18]:
{'old_min': 264.148,
'old_max': 520.245,
'old_mean': 347.9752,
'old_median': 350.69598,
'old_std': 49.78903,
'new_min': 264.706,
'new_max': 519.89703,
'new_mean': 356.3459,
'new_median': 364.74097,
'new_std': 53.15585,
'diff_min': -6.8740234,
'diff_max': 43.521973,
'diff_mean': 8.370939,
'diff_median': 4.911972,
'diff_std': 9.06017,
'roughness_min': 0.024993896484375,
'roughness_max': 16.99798583984375,
'roughness_mean': 0.8664416971654487,
'roughness_median': 0.501007080078125,
'roughness_std': 0.9008248692669407,
'slope_min': 0.00030909907704678905,
'slope_max': 83.09449329265449,
'slope_mean': 15.825287450127606,
'slope_median': 9.967977228482425,
'slope_std': 14.601294984565905,
'aspect_min': 0.0,
'aspect_max': 359.9997588896344,
'aspect_mean': 172.74772778897523,
'aspect_median': 183.00927436545985,
'aspect_std': 104.09906838438863,
'curvature_min': -1643.609619140625,
'curvature_max': 2639.605712890625,
'curvature_mean': 1.7315227064581844,
'curvature_median': 0.8087158203125,
'curvature_std': 43.72569863077058,
'terrain_ruggedness_index_min': 0.024669455802839076,
'terrain_ruggedness_index_max': 27.121265279417575,
'terrain_ruggedness_index_mean': 0.8939862986187255,
'terrain_ruggedness_index_median': 0.5406418612709043,
'terrain_ruggedness_index_std': 0.8938326265050937,
'rugosity_min': 1.000078699472661,
'rugosity_max': 9.8867222334156,
'rugosity_mean': 1.1083063933475281,
'rugosity_median': 1.0261556432873709,
'rugosity_std': 0.19958467669264315,
'profile_curvature_min': -1767.2495203822261,
'profile_curvature_max': 1433.1465250703304,
'profile_curvature_mean': -1.2548662332856713,
'profile_curvature_median': -0.5842965641336939,
'profile_curvature_std': 28.49170269374419,
'planform_curvature_min': -624.432616658446,
'planform_curvature_max': 978.0291066330788,
'planform_curvature_mean': 0.4766817029082818,
'planform_curvature_median': 0.18651537440086147,
'planform_curvature_std': 21.37279648106214}
In [19]:
REC = gpd.read_file("nzRec2_v5.gdb", engine='pyogrio', use_arrow=True)
REC
✔️ 3.4 s (2024-12-02T11:04:45/2024-12-02T11:04:48)
/home/ubuntu/.local/lib/python3.10/site-packages/pyogrio/raw.py:287: UserWarning: More than one layer found in 'nzRec2_v5.gdb': 'riverlines' (default), 'Hydro_Net_Junctions', 'rec2ws', 'rec2_rain_runoff_V5', 'rec2_stew_baserock_V5', 'rec2_stew_toprock_V5', 'rec2_utility_variables_V5', 'rec2_lcdb2_V5', 'rec2_lcdb1_V5', 'rec2_lcdb3_V5', 'rec2_ni_baserock_V5', 'rec2_ni_toprock_V5', 'rec2_si_baserock_V5', 'rec2_si_toprock_V5', 'rec2_elevband_rain_V5', 'LAYERKEYTABLE', 'APUNIQUEID', 'riverlines_FS', 'N_1_Desc', 'N_1_EDesc', 'N_1_EStatus', 'N_1_ETopo', 'N_1_FloDir', 'N_1_JDesc', 'N_1_JStatus', 'N_1_JTopo', 'N_1_JTopo2', 'N_1_Props'. Specify layer parameter to avoid this warning. with open_arrow(
Out[19]:
| OBJECTID | HydroID | NextDownID | CATAREA | CUM_AREA | nzsegment | Enabled | LENGTHDOWN | Headwater | Hydseq | StreamOrde | euclid_dis | upElev | downElev | upcoordX | downcoordX | downcoordY | upcoordY | sinuosity | nzreach_re | headw_dist | segslpmean | LID | reachtype | FROM_NODE | TO_NODE | Shape_Leng | Shape_Length | FLOWDIR | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 9 | 218553.786394 | 2.185538e+05 | 1000005 | 1 | 1801.260253 | 1 | 1 | 1 | 199.150697 | 119.994514 | 90.628487 | 1.601529e+06 | 1.601574e+06 | 6.192386e+06 | 6.192580e+06 | 1.072393 | 1000004 | 0 | 8.388180 | 0 | 0 | 1 | 2 | 213.567866 | 213.567866 | 1 | MULTILINESTRING ((1601528.857 6192580.4, 16015... |
| 1 | 2 | 2 | 9 | 455995.848576 | 4.559958e+05 | 1000003 | 1 | 1801.260253 | 1 | 2 | 1 | 537.313689 | 163.071075 | 90.628487 | 1.601783e+06 | 1.601574e+06 | 6.192386e+06 | 6.192881e+06 | 1.082775 | 1000005 | 0 | 7.678527 | 0 | 0 | 3 | 2 | 581.789877 | 581.789877 | 1 | MULTILINESTRING ((1601782.856 6192881.076, 160... |
| 2 | 3 | 3 | 4 | 461385.667778 | 4.613857e+05 | 1000008 | 1 | 335.656220 | 1 | 4 | 1 | 596.782205 | 64.632355 | 20.004650 | 1.599355e+06 | 1.599237e+06 | 6.191779e+06 | 6.192364e+06 | 1.163688 | 1000009 | 0 | 4.276652 | 0 | 0 | 4 | 5 | 694.468045 | 694.468045 | 1 | MULTILINESTRING ((1599355.244 6192363.83, 1599... |
| 3 | 4 | 4 | -1 | 135806.633398 | 2.330307e+06 | 1000010 | 1 | 0.000000 | 0 | 9 | 2 | 288.064229 | 20.004650 | 4.172440 | 1.599237e+06 | 1.598982e+06 | 6.191913e+06 | 6.191779e+06 | 1.165213 | 1000010 | 788 | 3.145852 | 0 | 0 | 5 | 6 | 335.656220 | 335.656220 | 1 | MULTILINESTRING ((1599237.074 6191778.666, 159... |
| 4 | 5 | 5 | -1 | 863421.445984 | 8.634214e+05 | 1000006 | 1 | 0.000000 | 1 | 5 | 1 | 1020.122542 | 143.468826 | 4.437851 | 1.602684e+06 | 1.602267e+06 | 6.191623e+06 | 6.192554e+06 | 1.182202 | 1000006 | 0 | 7.760943 | 0 | 0 | 7 | 8 | 1205.990464 | 1205.990464 | 1 | MULTILINESTRING ((1602683.557 6192553.928, 160... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 593512 | 593614 | 593513 | -1 | 51990.267885 | 1.174888e+09 | 4170001 | 1 | NaN | 0 | 124017 | 6 | 0.000000 | 0.000000 | 0.000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 85116 | 604841 | 542.278140 | 542.278140 | 1 | MULTILINESTRING ((1975828.498 5785372.292, 197... |
| 593513 | 593616 | 593514 | -1 | 150666.238187 | 2.061735e+08 | 1040001 | 1 | NaN | 0 | 27152 | 5 | 0.000000 | 12.272364 | 0.000000 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 26864 | 604842 | 1675.015974 | 1675.015974 | 1 | MULTILINESTRING ((1731848.082 6017486.385, 173... |
| 593514 | 593620 | 593515 | -1 | 53587.929573 | 5.503258e+07 | 7260002 | 1 | NaN | 0 | 239462 | 5 | 0.000000 | 4.917045 | 3.455193 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 244056 | 604843 | 156.868127 | 156.868127 | 1 | MULTILINESTRING ((1789499.342 5528919.92, 1789... |
| 593515 | 593621 | 593516 | 593515 | 8667.469549 | 5.432822e+07 | 7260001 | 1 | NaN | 0 | 239461 | 5 | 0.000000 | 6.187159 | 4.917045 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000 | 0 | 0 | 0.000000 | 0 | 0 | 244073 | 244056 | 116.148688 | 116.148688 | 1 | MULTILINESTRING ((1789607.824 5528878.42, 1789... |
| 593516 | 593622 | 593517 | -1 | 492400.109562 | 1.856727e+08 | 7247115 | 1 | 529.317107 | 0 | 249508 | 5 | 964.116694 | 5.069054 | 1.886685 | 1.782739e+06 | 1.782650e+06 | 5.496834e+06 | 5.497794e+06 | 1.086063 | 7047051 | 36796 | 0.272352 | 0 | 0 | 252501 | 604844 | 151.702564 | 151.702564 | 1 | MULTILINESTRING ((1782739.481 5497793.71, 1782... |
593517 rows × 30 columns
In [20]:
LCDB = gpd.read_file("lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gpkg", engine='pyogrio', use_arrow=True)
LCDB
✔️ 6.03 s (2024-12-02T11:04:48/2024-12-02T11:04:54)
Out[20]:
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | Wetland_18 | Wetland_12 | Wetland_08 | Wetland_01 | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb2000096676 | MULTIPOLYGON (((1613722.435 5425797.372, 16137... |
| 1 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Regional Council | 2019-12-01 | lcdb1000513359 | MULTIPOLYGON (((1816770.219 5947804.627, 18167... |
| 2 | Mangrove | Mangrove | Mangrove | Mangrove | Mangrove | 70 | 70 | 70 | 70 | 70 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb1000182160 | MULTIPOLYGON (((1715672.186 5958842.706, 17156... |
| 3 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | yes | yes | yes | yes | yes | no | no | no | no | no | Terralink | 2004-06-30 | lcdb1000065930 | MULTIPOLYGON (((1705330.918 6088979.74, 170531... |
| 4 | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | Estuarine Open Water | 22 | 22 | 22 | 22 | 22 | no | no | no | no | no | no | no | no | no | no | Regional Council | 2019-12-01 | lcdb1000065472 | MULTIPOLYGON (((1761684.636 5789742.527, 17616... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 511099 | Low Producing Grassland | Low Producing Grassland | Low Producing Grassland | Manuka and/or Kanuka | Low Producing Grassland | 41 | 41 | 41 | 52 | 41 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb1000505810 | MULTIPOLYGON (((1785112.194 5684560.595, 17851... |
| 511100 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb2000193885 | MULTIPOLYGON (((1607510.75 5432591.699, 160754... |
| 511101 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | 71 | 71 | 71 | 71 | 71 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb2000219027 | MULTIPOLYGON (((1603592.31 5269947.382, 160360... |
| 511102 | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | Herbaceous Freshwater Vegetation | 45 | 45 | 45 | 45 | 45 | yes | yes | yes | yes | yes | yes | yes | yes | yes | yes | Regional Council | 2014-06-30 | lcdb1000417326 | MULTIPOLYGON (((1822629.596 5477414.336, 18226... |
| 511103 | High Producing Exotic Grassland | High Producing Exotic Grassland | Exotic Forest | Exotic Forest | Exotic Forest | 40 | 40 | 71 | 71 | 71 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30 | lcdb1000145507 | MULTIPOLYGON (((1704402.741 5624819.729, 17044... |
511104 rows × 24 columns
In [21]:
LCDB[LCDB.intersects(largest_area.geometry)]
✔️ 564 ms (2024-12-02T11:04:55/2024-12-02T11:04:55)
Out[21]:
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | Wetland_18 | Wetland_12 | Wetland_08 | Wetland_01 | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 17220 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | High Producing Exotic Grassland | 71 | 71 | 71 | 71 | 40 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000053082 | MULTIPOLYGON (((2054973.343 5805787.881, 20549... |
| 18662 | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | Gravel or Rock | 16 | 16 | 16 | 16 | 16 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30 | lcdb1000016122 | MULTIPOLYGON (((2048258.705 5805949.487, 20482... |
| 20295 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127897 | MULTIPOLYGON (((2054496.624 5803468.778, 20545... |
| 23069 | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | Broadleaved Indigenous Hardwoods | 54 | 54 | 54 | 54 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb1000119141 | MULTIPOLYGON (((2054140.516 5804625.89, 205416... |
| 349500 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127894 | MULTIPOLYGON (((2056301.039 5802553.702, 20562... |
| 388204 | Exotic Forest | Exotic Forest | Exotic Forest | Exotic Forest | Broadleaved Indigenous Hardwoods | 71 | 71 | 71 | 71 | 54 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2011-06-30 | lcdb1000127900 | MULTIPOLYGON (((2054727.599 5803954.788, 20547... |
| 445455 | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Mixed Exotic Shrubland | Gravel or Rock | 56 | 56 | 56 | 56 | 16 | no | no | no | no | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01 | lcdb2000515912 | MULTIPOLYGON (((2055082.956 5803253.022, 20550... |
In [30]:
LCDB[LCDB.intersects(largest_area.geometry)].drop(columns="EditDate").explore("Name_2018", legend=True, tiles="ESRI.WorldImagery")
✔️ 368 ms (2024-12-02T11:10:54/2024-12-02T11:10:54)
Out[30]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [23]:
LCDB.loc[LCDB.intersects(largest_area.geometry), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True})
✔️ 290 ms (2024-12-02T11:04:58/2024-12-02T11:04:58)
Out[23]:
Class_2018 71 Wetland_18 False Onshore_18 True Name: 0, dtype: object
In [24]:
all_streams = REC.unary_union
✔️ 10.4 s (2024-12-02T11:04:58/2024-12-02T11:05:08)
/tmp/ipykernel_3591039/3529544978.py:1: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead. all_streams = REC.unary_union
In [25]:
features["distance_to_river"] = largest_area.geometry.distance(all_streams)
features["distance_to_river"]
✔️ 156 ms (2024-12-02T11:05:09/2024-12-02T11:05:09)
Out[25]:
0.0
In [26]:
def get_features(row):
geom = row.geometry
masked_old, transform = mask(old, [geom], nodata=np.nan)
nonan = np.argwhere(~np.isnan(masked_old[0]))
top_left = nonan.min(axis=0)
bottom_right = nonan.max(axis=0)
masked_old = masked_old[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
masked_new, transform = mask(new, [geom], nodata=np.nan)
masked_new = masked_new[0][top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
raster = rasterize([geom], out_shape=result.shape, transform=new.transform)
masked_diff = np.where(raster, diff, np.nan)
masked_diff = masked_diff[top_left[0]:bottom_right[0]+1, top_left[1]:bottom_right[1]+1]
row = row.to_dict()
row.update({
"old_min": np.nanmin(masked_old),
"old_max": np.nanmax(masked_old),
"old_mean": np.nanmean(masked_old),
"old_median": np.nanmedian(masked_old),
"old_std": np.nanstd(masked_old),
"new_min": np.nanmin(masked_new),
"new_max": np.nanmax(masked_new),
"new_mean": np.nanmean(masked_new),
"new_median": np.nanmedian(masked_new),
"new_std": np.nanstd(masked_new),
"diff_min": np.nanmin(masked_diff),
"diff_max": np.nanmax(masked_diff),
"diff_mean": np.nanmean(masked_diff),
"diff_median": np.nanmedian(masked_diff),
"diff_std": np.nanstd(masked_diff),
"distance_to_river": geom.distance(all_streams)
})
row.update(LCDB.loc[LCDB.intersects(geom), ["Class_2018", "Wetland_18", "Onshore_18"]].mode().iloc[0].replace({"no": False, "yes": True}))
attribute_names = ["roughness", "slope", "aspect", "curvature", "terrain_ruggedness_index", "rugosity", "profile_curvature", "planform_curvature"]
attributes = xdem.terrain.get_terrain_attribute(
masked_diff,
resolution=old.res,
attribute=attribute_names
)
for i, name in enumerate(attribute_names):
row.update({
f"{name}_min": np.nanmin(attributes[i]),
f"{name}_max": np.nanmax(attributes[i]),
f"{name}_mean": np.nanmean(attributes[i]),
f"{name}_median": np.nanmedian(attributes[i]),
f"{name}_std": np.nanstd(attributes[i]),
})
return pd.Series(row)
#get_features(areas.iloc[0])
features = areas.progress_apply(get_features, axis=1)
features
✔️ 5 min 39 s (2024-12-02T11:05:09/2024-12-02T11:10:49)
0%| | 0/188 [00:00<?, ?it/s]
Out[26]:
| geometry | area | old_min | old_max | old_mean | old_median | old_std | new_min | new_max | new_mean | new_median | new_std | diff_min | diff_max | diff_mean | diff_median | diff_std | distance_to_river | Class_2018 | Wetland_18 | Onshore_18 | roughness_min | roughness_max | roughness_mean | roughness_median | roughness_std | slope_min | slope_max | slope_mean | slope_median | slope_std | aspect_min | aspect_max | aspect_mean | aspect_median | aspect_std | curvature_min | curvature_max | curvature_mean | curvature_median | curvature_std | terrain_ruggedness_index_min | terrain_ruggedness_index_max | terrain_ruggedness_index_mean | terrain_ruggedness_index_median | terrain_ruggedness_index_std | rugosity_min | rugosity_max | rugosity_mean | rugosity_median | rugosity_std | profile_curvature_min | profile_curvature_max | profile_curvature_mean | profile_curvature_median | profile_curvature_std | planform_curvature_min | planform_curvature_max | planform_curvature_mean | planform_curvature_median | planform_curvature_std | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 176 | POLYGON ((2054419 5803910, 2054419 5803909, 20... | 530308.0 | 264.148010 | 520.244995 | 347.975189 | 350.695984 | 49.789028 | 264.705994 | 519.897034 | 356.345886 | 364.740967 | 53.155849 | -6.874023 | 43.521973 | 8.370939 | 4.911972 | 9.060170 | 0.000000 | 71 | False | True | 0.024994 | 16.997986 | 0.866442 | 0.501007 | 0.900825 | 0.000309 | 83.094493 | 15.825287 | 9.967977 | 14.601295 | 0.000000 | 359.999759 | 172.747728 | 183.009274 | 104.099068 | -1643.609619 | 2639.605713 | 1.731523 | 0.808716 | 43.725699 | 0.024669 | 27.121265 | 0.893986 | 0.540642 | 0.893833 | 1.000079 | 9.886722 | 1.108306 | 1.026156 | 0.199585 | -1767.249520 | 1433.146525 | -1.254866 | -0.584297 | 28.491703 | -624.432617 | 978.029107 | 0.476682 | 0.186515 | 21.372796 |
| 155 | POLYGON ((2054514 5804756, 2054514 5804755, 20... | 319463.0 | 384.855011 | 949.697021 | 698.896851 | 715.611023 | 132.071274 | 384.348999 | 949.190002 | 684.987671 | 692.672974 | 129.693771 | -52.056030 | 5.609009 | -13.909487 | -9.703979 | 12.755947 | 0.000000 | 54 | False | True | 0.027039 | 10.799988 | 1.465731 | 1.246033 | 0.944709 | 0.021028 | 77.297658 | 26.350639 | 24.870094 | 14.383982 | 0.000000 | 359.999410 | 183.019811 | 185.773786 | 100.256831 | -954.620361 | 692.401123 | -1.137023 | -1.690674 | 58.625323 | 0.035931 | 10.752627 | 1.490565 | 1.282474 | 0.923966 | 1.000176 | 4.733187 | 1.209573 | 1.132283 | 0.229804 | -440.434936 | 585.141157 | 0.664059 | 0.707163 | 35.661075 | -369.479204 | 312.140260 | -0.473099 | -0.604222 | 31.241830 |
| 37 | POLYGON ((2056589 5808659, 2056589 5808657, 20... | 304243.0 | 543.239014 | 968.459961 | 752.026978 | 758.137024 | 89.943939 | 542.731018 | 968.713989 | 750.921082 | 757.086975 | 90.210747 | -13.114990 | 5.141968 | -1.105809 | -0.944031 | 1.355662 | 0.000000 | 12 | False | True | 0.014038 | 6.578003 | 0.709351 | 0.576965 | 0.492890 | 0.019821 | 66.326472 | 13.694865 | 11.374013 | 9.495378 | 0.000000 | 359.998682 | 180.263642 | 180.474840 | 104.744639 | -472.998047 | 416.595459 | -1.971487 | -1.379395 | 39.778407 | 0.014390 | 5.968466 | 0.742228 | 0.609774 | 0.502908 | 1.000017 | 2.568007 | 1.063572 | 1.033145 | 0.087242 | -275.317190 | 316.405986 | 1.262352 | 0.660475 | 23.922942 | -335.581215 | 226.580603 | -0.709152 | -0.401965 | 21.456914 |
| 115 | POLYGON ((2054495 5806342, 2054495 5806340, 20... | 173781.0 | 333.551025 | 461.656006 | 378.365387 | 375.186005 | 26.501810 | 334.056000 | 461.194000 | 381.083679 | 378.898010 | 26.786768 | -4.665009 | 12.177979 | 2.718257 | 2.453003 | 1.906919 | 0.000000 | 12 | False | True | 0.028015 | 6.747986 | 0.664326 | 0.454987 | 0.605888 | 0.028203 | 69.110675 | 12.514475 | 8.870377 | 11.030167 | 0.000000 | 359.999000 | 164.568461 | 177.295377 | 100.462405 | -553.802490 | 687.600708 | 2.607275 | 1.797485 | 41.001715 | 0.032552 | 7.412247 | 0.699980 | 0.489235 | 0.623502 | 1.000124 | 2.874868 | 1.065598 | 1.021682 | 0.119312 | -479.441262 | 366.685717 | -2.050307 | -1.421865 | 27.232389 | -357.169244 | 403.752029 | 0.556968 | 0.261994 | 19.498026 |
| 126 | POLYGON ((2054821 5806083, 2054821 5806082, 20... | 147413.0 | 405.266998 | 984.378967 | 666.230042 | 680.250000 | 128.762848 | 404.644012 | 984.406982 | 662.956543 | 677.700012 | 128.024765 | -25.545959 | 10.151001 | -3.273446 | -2.276978 | 3.636282 | 0.000000 | 71 | False | True | 0.049988 | 8.898010 | 1.275577 | 1.083984 | 0.813165 | 0.040820 | 73.866107 | 23.143639 | 21.255515 | 13.418823 | 0.000000 | 359.998574 | 184.088882 | 176.004473 | 102.389142 | -767.303467 | 707.580566 | -4.001010 | -3.497314 | 62.809494 | 0.051240 | 8.480720 | 1.322321 | 1.132742 | 0.821058 | 1.000402 | 3.628650 | 1.171587 | 1.106380 | 0.191066 | -536.233840 | 418.001613 | 2.441112 | 1.776779 | 38.928714 | -393.931141 | 379.570854 | -1.559898 | -1.174689 | 33.025140 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 51 | POLYGON ((2057180 5807155, 2057180 5807154, 20... | 4309.0 | 754.164001 | 862.838013 | 810.229065 | 812.798035 | 27.570267 | 753.013000 | 862.265991 | 808.988464 | 811.622009 | 27.611164 | -5.234009 | 0.836975 | -1.240615 | -1.079041 | 0.706411 | 377.665475 | 69 | False | True | 0.091980 | 3.499023 | 0.810893 | 0.724518 | 0.466198 | 0.258646 | 54.086741 | 15.535889 | 14.296695 | 9.054989 | 0.000000 | 359.873676 | 180.929910 | 181.396905 | 105.351676 | -295.800781 | 442.993164 | -7.090177 | -4.492188 | 47.193806 | 0.086022 | 3.618491 | 0.855655 | 0.758100 | 0.486536 | 1.001061 | 1.724303 | 1.077381 | 1.049964 | 0.088696 | -238.362703 | 175.115296 | 2.803598 | 1.571423 | 27.481605 | -176.705037 | 204.630461 | -4.286579 | -2.724701 | 26.676936 |
| 77 | POLYGON ((2057221 5806249, 2057221 5806248, 20... | 4231.0 | 463.237000 | 558.276001 | 506.580353 | 501.701996 | 21.162735 | 462.401001 | 557.640991 | 504.420471 | 500.536987 | 21.546041 | -9.316010 | 0.420044 | -2.159984 | -1.524994 | 1.754438 | 0.000000 | 16 | False | True | 0.123962 | 5.431030 | 1.354514 | 1.208984 | 0.781550 | 0.243961 | 63.330647 | 24.899263 | 23.870917 | 12.951836 | 0.117864 | 359.932983 | 184.918589 | 175.719571 | 102.218871 | -365.005493 | 229.119873 | -11.565558 | -9.808350 | 54.859841 | 0.102363 | 4.997128 | 1.382453 | 1.254762 | 0.757155 | 1.000989 | 2.266732 | 1.179069 | 1.126430 | 0.172435 | -180.493919 | 211.376430 | 5.078110 | 3.214066 | 31.870435 | -207.633207 | 137.569661 | -6.487448 | -5.183216 | 31.373435 |
| 181 | POLYGON ((2054921 5802132, 2054921 5802131, 20... | 4168.0 | 727.848999 | 777.617981 | 750.956482 | 751.011475 | 10.164408 | 727.346008 | 777.073975 | 748.487061 | 747.432495 | 10.804439 | -7.746033 | 0.098999 | -2.469438 | -1.833984 | 1.776451 | 313.839735 | 71 | False | True | 0.050049 | 2.488037 | 0.713715 | 0.641541 | 0.372780 | 0.161113 | 41.709637 | 14.652946 | 13.619203 | 7.531259 | 0.239895 | 359.927197 | 183.505383 | 172.413544 | 104.640565 | -209.301758 | 189.300537 | -3.036568 | -1.809692 | 25.295615 | 0.050731 | 2.201049 | 0.710312 | 0.651643 | 0.352425 | 1.000237 | 1.344907 | 1.051155 | 1.036017 | 0.048002 | -116.756886 | 137.333515 | 1.427443 | 0.626500 | 14.888978 | -90.982558 | 93.936800 | -1.609124 | -0.725124 | 13.837087 |
| 112 | POLYGON ((2056165 5805098, 2056165 5805097, 20... | 4168.0 | 330.009003 | 414.976013 | 347.387390 | 338.135986 | 19.942204 | 329.310974 | 414.473022 | 345.676178 | 335.724487 | 20.337126 | -6.242035 | 0.404999 | -1.711227 | -1.103989 | 1.331259 | 20.061887 | 16 | False | True | 0.067993 | 5.648010 | 0.916036 | 0.627975 | 0.812676 | 0.072869 | 67.422201 | 16.667242 | 11.679819 | 14.067467 | 0.000371 | 359.989076 | 160.836444 | 161.757142 | 101.798624 | -645.401001 | 310.797119 | -10.793974 | -5.749512 | 54.868025 | 0.061296 | 7.057928 | 0.977414 | 0.664636 | 0.868992 | 1.000381 | 2.796020 | 1.118655 | 1.040430 | 0.184431 | -213.102706 | 411.064214 | 7.399746 | 3.237263 | 38.168409 | -249.806248 | 200.170430 | -3.394228 | -2.160090 | 26.729371 |
| 70 | POLYGON ((2056954 5806554, 2056954 5806551, 20... | 4069.0 | 580.815002 | 632.517029 | 607.255920 | 607.020020 | 11.594285 | 580.307007 | 632.016968 | 605.696106 | 604.859009 | 11.600400 | -5.670044 | -0.108032 | -1.559806 | -1.331970 | 0.911266 | 44.852111 | 69 | False | True | 0.078003 | 2.404053 | 0.727969 | 0.664001 | 0.367533 | 0.283751 | 40.625595 | 14.404642 | 13.091004 | 7.748489 | 0.185128 | 359.981782 | 184.885313 | 190.514454 | 104.413637 | -207.397461 | 179.193115 | -6.322378 | -4.388428 | 34.206604 | 0.078687 | 2.184222 | 0.748892 | 0.687885 | 0.359458 | 1.000577 | 1.326722 | 1.056693 | 1.041069 | 0.049973 | -87.630775 | 112.695939 | 2.519276 | 1.309808 | 19.584720 | -115.848622 | 91.562340 | -3.803102 | -2.074324 | 19.314333 |
188 rows × 61 columns
In [27]:
features.crs = new.crs
✔️ 157 ms (2024-12-02T11:10:49/2024-12-02T11:10:50)
In [28]:
features.explore("Class_2018", legend=True, tiles="ESRI.WorldImagery")
✔️ 1.95 s (2024-12-02T11:10:50/2024-12-02T11:10:52)
Out[28]:
Make this Notebook Trusted to load map: File -> Trust Notebook